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Abstract 

This article gives a direct formula for the computation of B (n) using the asymptotic 
formula 



where n is even and n ^ 1. This is simply based on the fact that C (n) is very near 
1 when n is large and since B (n) = 2 ^^S^tt exactly. The formula chosen for the Zeta 
function is the one with prime numbers from the well-known Euler product for (n). 
This algorithm is far better than the recurrence formula for the Bernoulli numbers even 
if each B{n) is computed individually. The author could compute B (750, 000) in a few 
hours. The current record of computation is now (as of Feb. 2007) B (5, 000, 000) a 
number of (the numerator) of 27332507 decimal digits is also based on that idea. 



1 The need for a single computation 

This algorithm came once in 1996 when the authors wanted to compute large Bernoulli num- 
bers using a well-known computer algebra system system like Maple or Mathematica. These 
programs used Faulhaber's recurrence [2, 5] formula which is nice but unsuitable for large 
computations. We quickly came to the conclusion that 5(10000) was out of reach even with 
a powerful computer. This is where we realized that for n large the actual formula is simply 
B (n) = 2 where n is even and not counting the sign, for n=1000 the approximation is 
good to more than 300 decimal digits where B (1000) is of the order of 1770 digits. To carry 
out the exact computation of B (1000) one has only to compute first the principal term in the 
asymptotic formula and secondly just a few terms in the Euler product (up to p = 59). The 
second idea was that the fractional part of the Bernoulli numbers can also be computed very 
fast with the help of the Von Staudt-Clausen formula. So finally, the need is only to compute 
Bn with enough precision so that the remainder is < 1 and apply the Von Staudt-Clausen 
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formula for the fractional part to finally add the 2 results. Note : Mathematica now uses a 
much more efficient algorithm partly due to these results presented here. 



2 The Von Staudt- Clausen formula 



The formula is, for k > 1, 



mod 1. 



The sum being extended over primes p such that {p — l)\2k [5]. In other words, for -B(IO) 
the sum is 



In terms of computation, when n is of the order of 1000000 it goes very fast to compute the 
fractional part of Bn- The only thing that remains to be done then is the principal part or 
integer part of Bn- 

3 The Euler product 

The Euler product of the zeta function is 



Where s > 1 and p is prime. This is the error term in For any given n there are 
primes compared to n. Translated into the program it means less operations to carry, the 
program stops when p'^ is of the order of B (n). 



4 The final program 

The Maple program uses a high precision value of 2tt and a routine for the Von Staudt- 
Clausen formula. That program held the record of the computation of Bernoulli Numbers 
from 1996 to 2002, after that others made more efficient programs using C++ and high 
precision packages like Kellner and Pavlyk (see table 1) and could reach 5(5,000,000). 
The program was used in 2003 to verify Agoh's conjecture up to n=49999 by the authors. 
Agoh's conjecture is 



5(10) = 1 - 1/2 - 1/3 - 1/11 = 5/66. 




1 mod p 
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is true iff p is prime. The congruence is not obvious since p-Bp-i is a fraction. Tlie standard 
method reduces first the numerator mod p, then re-evaluates the fraction, then reduces the 
numerator mod p. The final fraction is always smaller than 1 and the result of a/b mod p 
is solved by finding k such that a = bk mod p. There are 3 parts in the main program 
which may take time. First the computation of (27r)" and n\. Secondly, the evaluation of 
the Von Staudt-Clausen formula and thirdly the computation of the Euler product. On a 
medium sized computer (Pentium 2.4 Ghz with Maple 10 and 1 gigabyte of memory). The 
run time for 5(20000) is about 9 seconds and the number is 61382 digits long including 1 sec- 
ond to read the value of vr to high-precision from the disk. Here are the timings for that run : 

- Product with primes up to 1181 at 61382 digits of precision : 7 seconds. 

- Exponentiation of 27r and n\ : less than 1 second. 

- Computation of 20000! : negligible. 

- Computation of Von Staudt-Clausen expression : negligible. 

When n increases the time taken to evaluate the product with primes is what takes the most. 
A value of vr to several thousands digits is necessary. Maple can supply many thousands but 
a file containing 1 million is easily found on the internet and is much faster. In this program 
TT is renamed pi with no capitals. The Bernoulli numbers up to n = 100 are within the 
program mainly for speed when n is small. 



BERN:=proc(n: : integer) 

local d, z, oz, i, p, pn, pnl, f, s, pi, tl, t2; 
global Digits; 

Iprint (' start at time' = timeO); 

if n = 1 then -1/2 

elif n = then 1 

elif n < then ERROR (' argument must be >= 0') 
elif irem(n, 2) = 1 then 

elif n <= 100 then op(iquo(n, 2), [1/6, -1/30, 1/42, -1/30, 
5/66, -691/2730, 7/6, -3617/510, 43867/798, -174611/330, 
854513/138, -236364091/2730, 8553103/6, -23749461029/870, 
8615841276005/14322, -7709321041217/510, 2577687858367/6, 
-26315271553053477373/1919190, 2929993913841559/6, 
-261082718496449122051/13530 , 1520097643918070802691/1806 , 
-27833269579301024235023/690, 596451111593912163277961/282, 
-5609403368997817686249127547/46410, 
495057205241079648212477525/66 , 
-801165718135489957347924991853/1590, 
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29149963634884862421418123812691/798, 

-2479392929313226753685415739663229/870, 

84483613348880041862046775994036021/354, 

-1215233140483755572040304994079820246041491/56786730, 

12300585434086858541953039857403386151/6, 

-106783830147866529886385444979142647942017/510, 

1472600022126335654051619428551932342241899101/64722, 

-78773130858718728141909149208474606244347001/30, 

1505381347333367003803076567377857208511438160235/4686, 

-5827954961669944110438277244641067365282488301844260429/ 

140100870, 

34152417289221168014330073731472635186688307783087/6, 
-24655088825935372707687196040585199904365267828865801/30, 
414846365575400828295179035549542073492199375372400483487/ 
3318 , -46037842994794576469355749690190468497942578727512\ 
88919656867/230010 , 1677014149185145836823154509786269900\ 
207736027570253414881613/498, -20245761959352903602311311\ 
60111731009989917391198090877281083932477/3404310, 660714\ 
61941767865357384784742626149627783068665338893176199698\ 
3/6, -131142648867401750799551142401931184334575027557202\ 
8644296919890574047/61410 , 117905727902108279988412335124\ 
9215083775254949669647116231545215727922535/272118, -1295\ 
58594820753752798942782853857674965934148371943514302331\ 
6326829946247/1410, 1220813806579744469607301679413201203\ 
958508415202696621436215105284649447/6 , -2116004495972665\ 
13097597728109824233673043954389060234150638733420050668\ 
349987259/4501770 , 67908260672905495624051 1 17546403605607\ 
342195728504487509073961249992947058239/6 , -9459803781912\ 
21252952274330694937218727028415330669361333856962043113\ 
95415197247711/33330] ) 

d := 4 

+ truncCevalhf ((InGAMMACn + 1) - n*ln(2*Pi) )/ln(10) ) ) 

+ length (n) ; 
Iprint ( 'using ' . d . ' Digits'); 
s := truncCevalhf (exp(0.5*d*ln(10)/n))) + 1; 
Digits := d; 

P := 1; 

tl := 1. ; 
t2 := tl; 



4 



Iprint (' start small prime loop at time' = timeO); 
while p <= s do 

p := nextprime(p) ; 



pn 


:= p n; 


pnl 


:= pn - 


tl 


:= pn*tl; 


t2 


:= pnl*t2 


end do; 




gcO; 





Iprint (status) ; 

Iprint ('used primes up to and including ' . p) ; 
Iprint ( 'finish small prime loop at time' = timeO); 
z := tl/t2; 
gcO; 

Iprint (status) ; 

Iprint ( 'finish full prec. division at time' = timeO); 
oz : = ; 

while oz <> z do 
oz := z; 

p := nextprime(p) ; 

Digits := max(d - iloglO(pn), 9); 

pn := Float (p,0) ; 

pn := p"n; 

pnl := z/pn; 

Digits := d; 

z := z + pnl 
end do; 
gcO; 

Iprint (status) ; 

Iprint ('used primes up to and including ' . p) ; 
Iprint ( 'finish big prime loop at time' = timeO); 
p : = evalf (2*pi) ; 
gcO; 

Iprint (status) ; 

Iprint ( 'finish 2*Pi at time' = timeO); 
f := n!; 

gcO; 

Iprint (status) ; 

Iprint ( 'finish factorial at time' = timeO); 
pn := p~n; 
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gcO; 

Iprint (status) ; 

Iprint ( 'finish (2*Pi)"n at time' = timeO); 

z := 2*z*f/pn; 

gcO; 

Iprint (status) ; 
Iprint ( 

'finish 2*z*n ! / (2*Pi) "n (multiply and divide) at time' 
= timeO) ; 
s := 0; 

for p in numtheory [divisors] (n) do 

pi := p + 1; if isprime(pl) then s := s + 1/pl end if 
end do; 
gcO; 

Iprint (status) ; 

Iprint ( 'finish divisors of n loop at time' = timeO); 

s := frac(s) ; 

if irem(n, 4) = then 

if s < 1/2 then z := -round(z) - s 

else z := -trunc(z) - s 

end if 

else 

s := 1 - s; 

if s < 1/2 then z := round (z) + s 

else z := trunc(z) + s 

end if 
end if ; 
gcO; 

Iprint (status) ; 

Iprint ('done at time' = timeO); 
z 

end if 

end: 
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Table 1: History of the computation of Bernoulli numbers 
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